 % MATLAB 代码：基于 SCNR 和 SINR 的优化问题

% 参数初始化
P_T = 1;                   % 总发射功率
alpha = 0.5;                % 用于权衡感知和通信的权重系数
sigma_s = 0.1;              % 感知噪声功率
sigma_c = 0.1;              % 通信噪声功率
gamma = 0.2;                % 隐蔽约束阈值

% 随机初始化信道向量
h_s = randn(2,1) + 1i*randn(2,1);  % 感知信道向量
h_c = randn(2,1) + 1i*randn(2,1);  % 通信信道向量

% 随机初始化杂波协方差矩阵
C = randn(2) + 1i*randn(2);
C = C'*C;  % 确保是半正定的

% 可行的 MA 位置约束（位置初始化）
X = [0; 0];          % 可移动天线的初始位置
X_feasible = [-1, 1; -1, 1];  % MA 位置的可行区域

% 迭代参数
max_iters = 100;     % 最大迭代次数
tolerance = 1e-3;    % 收敛容差

% 初始化波束成形向量和人工噪声向量
w = sqrt(P_T/2) * (randn(2,1) + 1i*randn(2,1));  % 随机初始化波束成形向量
z = sqrt(P_T/2) * (randn(2,1) + 1i*randn(2,1));  % 随机初始化人工噪声向量

% 初始化上一轮迭代的值以检查收敛
w_prev = w;
z_prev = z;
X_prev = X;

% 记录目标函数的值
objective_vals = [];

% 迭代凸逼近循环
for iter = 1:max_iters
    % 计算感知的 SCNR
    SCNR_s = (abs(h_s' * w)^2) / (w' * C * w + abs(h_s' * z)^2 + sigma_s^2);
    
    % 计算通信的 SINR
    SINR_c = (abs(h_c' * w)^2) / (abs(h_c' * z)^2 + sigma_c^2);
    
    % 目标函数：alpha * SCNR_s + (1 - alpha) * SINR_c
    objective = alpha * SCNR_s + (1 - alpha) * SINR_c;
    objective_vals = [objective_vals, objective];  % 记录目标函数值
    
    % 使用凸优化包 cvx 进行波束成形和人工噪声向量的优化
    cvx_begin quiet
        variable w_new(2) complex
        variable z_new(2) complex
        variable t_w
        variable t_z            
        % 优化目标函数
        maximize(alpha * (pow_pos(abs(h_s' * w), 2) / (pow_pos(abs(h_s' * z), 2) + sigma_s^2)) + ...
            (1 - alpha) * (pow_pos(abs(h_c' * w), 2) / (pow_pos(abs(h_c' * z), 2) + sigma_c^2)))
                  
        % 约束条件：总功率约束
        subject to
            t_w >= norm(w_new);
            t_z >= norm(z_new);
            t_w^2 + t_z^2 <= P_T;
%             norm(w_new)^2 + norm(z_new)^2 <= P_T
    cvx_end
    
    % 更新波束成形向量和人工噪声向量
    w = w_new;
    z = z_new;
    
    % 更新 MA 位置（使用投影梯度下降法）
    X_gradient = gradient_cov_constraint(h_s, h_c, w, z, X);  % 自定义梯度计算函数
    X = X - 0.1 * X_gradient;  % 梯度下降步骤
    
    % 将 MA 位置投影回可行区域
    X = max(min(X, X_feasible(:,2)), X_feasible(:,1));  % 将 X 限制在可行边界内
    
    % 收敛检查
    if norm(w - w_prev) < tolerance && norm(z - z_prev) < tolerance && norm(X - X_prev) < tolerance
        disp(['迭代在第 ', num2str(iter), ' 轮收敛']);
        break;
    end
    
    % 更新前一轮的值
    w_prev = w;
    z_prev = z;
    X_prev = X;
end

% 目标函数进展的绘图
figure;
plot(1:length(objective_vals), objective_vals);
xlabel('迭代次数');
ylabel('目标函数值');
title('目标函数进展');

% 输出结果
disp('最优波束成形向量 (w):');
disp(w);
disp('最优人工噪声向量 (z):');
disp(z);
disp('最优 MA 位置 (X):');
disp(X);

